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Abstract 

Two propulsion aeroelasticity codes were used to model the 
aeroelastic characteristics of an experimental forward-swept 
fan that encountered flutter during wind tunnel testing. Both of 
these three-dimensional codes model the unsteady flowfield 
due to blade vibrations using the Navier-Stokes equations. In 
the first approach, the unsteady flow equations are solved 
using an implicit time-marching approach. In the second 
approach, the unsteady flow equations are converted to a 
harmonic balance form and solved using a pseudo-time 
marching method. This paper describes the flutter calculations 
and compares the results to experimental measurements. 

Introduction 

Aircraft engine turbomachinery blades are susceptible to 
aeroelastic vibration problems and high-cycle fatigue failures. 
Flutter stability is likely to become a major challenge as 
aggressive new designs of fan, compressor and turbine blades 
are developed to reduce noise, improve performance, and 
reduce weight. Hence, it is important to develop and validate 
numerical tools that can be used to verily aeroelastic stability 
using high-fidelity physics-based models. Such numerical 
tools will enable gas turbine designers to develop new 
turbomachinery blading that will not flutter, thus improving 
safety, reducing development cycle time and cost, and 
enabling the targeted improvements listed earlier. 

Research has been on-going in the development, validation 
and application of high-fidelity models for aeroelastic 
vibrations in aircraft engine fan, compressor, and turbine 
blades (Ref. 1). Recent work has included time-domain 
solution of the Reynolds-averaged Navier-Stokes (RANS) 
equations to provide the unsteady flowfield and unsteady 
aerodynamic forces on the blades. An example of such work is 
the TURBO aeroelastic analysis code (Refs. 2 and 3). Such 
high-fidelity time-domain models require large numbers of 
computations and a long time for startup transients to decay 
before the final periodic solution is obtained. A second 


approach is to use the periodicity in time of typical 
turbomachinery flows to represent each flow variable by a 
Fourier series in time, leading to a harmonic balance form of 
the Navier-Stokes equations (Ref. 4). Solutions to these 
equations can be obtained using methods that are typically 
used for steady flow problems such as pseudo-time marching 
and local time stepping. Thus, the harmonic balance approach 
leads to a method that can be significantly faster than the 
typical time-domain solution method (Ref. 4). 

In the present study, the configuration selected is an 
experimental fan (Ref. 5) for which wind-tunnel 

measurements of performance and flutter vibrations are 
available. This fan was designed with aggressive goals for 
performance and noise reduction. During wind-tunnel testing, 
the fan performed well at design speed, and was successfully 
throttled to the stall line. However, flutter was encountered 
just above the operating line at part-speed conditions. The 
flutter mode was identified as the first bending mode of the 
airfoil, in a two nodal diameter forward-traveling wave 
pattern. In this paper, a comparison is presented between two 
flutter analyses applied to this fan: a time-domain RANS 
aeroelastic code (Ref. 2) (TURBO) and a harmonic balance 
(HB) RANS aeroelastic analysis and computer code developed 
by Hall et al. (Refs. 4 and 6). Both analyses solve the three- 
dimensional unsteady, Reynolds-Averaged Navier-Stokes 
equations with the ability to model a rotating blade row with 
harmonic blade vibrations or incoming periodic distortions. 
For flutter calculations, the blade vibration is prescribed to be 
the modal deflection and a frequency, both of which are 
calculated from a separate structural dynamics analysis. 
Computations are performed in a single blade passage for both 
steady and unsteady analyses. In the time-domain analysis 
(Ref. 2), the unsteady pressures are used to calculate a work- 
per-cycle that is then used to calculate the aerodynamic 
damping, which is used to determine aeroelastic stability in 
the selected vibration mode and traveling wave pattern. In the 
harmonic balance analysis, the resulting unsteady pressures 
are used to calculate a generalized aerodynamic force and then 
an eigenvalue problem is solved to calculate the aerodynamic 
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damping, which is used to determine stability of the blade in a 
specific vibration mode (Ref. 7). 

Note that both the aeroelastic codes were run with small 
amplitude vibrations of the fan blade to calculate linearized 
unsteady aerodynamics for a conventional linear flutter 
analysis. However, both the analysis procedures can also be 
used for larger amplitudes of vibration. 

Analysis 

Aeroelastic Model 

The equations of motion for a fan blade (with all blades 
assumed to be identical) can be written as 

[M\q] + [K\ q }=[A\ q } (1) 

where [M] and [K\ are generalized mass and stiffness 
matrices, {q} is the generalized displacement vector, and [A] 
is the blade vibration-dependent generalized aerodynamic 
force matrix. The matrices [M], [£], and [A] are of size 
NMxNM ; {q} is of size NMxl; NM is the number of modes. 

The elements of [M] and [£] are obtained from a free- 
vibration analysis using a commercial structural dynamics 
analysis software. The matrices [M] and [/£] are diagonal and 
their non-zero elements are related as 

K^M^ 2 (1 + 2%) ( 2 ) 

where co ? is the natural frequency of the I th mode, and £ is the 
structural damping ratio; usually the mode shapes are mass- 
normalized and therefore M ( = 1 . 

Since all the blades are identical (that is, a tuned rotor), the 
aeroelastic modes consist of individual blades vibrating with 
equal amplitudes at a fixed interblade phase angle between 
adjacent blades. Hence, the motion of the s th blade in t* h 
interblade phase angle mode can be written as 

(3) 

where co is the vibration frequency, <j r is the interblade phase 
angle related to nodal diameter (ND) pattern of the traveling 
wave and number of blades Abiades as 

a r = 2n ND / N blades (4) 

Thus, the equations of motion for a blade become 

-<* 2 [M]{q r }+[Khr} = [A r br} (5) 

The following subsections describe the flutter calculation 
method used with the harmonic balance code and separately 
with the time-domain TURBO code. 


HB Flutter Analysis 

The HB flutter analysis requires calculation of elements of 
the generalized aerodynamic force matrix [A r \. Unsteady 
flowfield computations are carried out for each vibration mode 
and assumed frequency. For a selected value of interblade 
phase angle, the harmonic balance code is used to calculate the 
unsteady pressure distribution on the blade surface, which is 
further used to calculate the (complex-valued) elements of the 
generalized aerodynamic force matrix [A r ]. This calculation is 
repeated for iV b i ades interblade phase angles given by 
Equation (4). 

To calculate flutter stability, Equation (5) is written in a 
standard eigenvalue form as: 

tt^HelkMo} (6) 

where 

[/ > ]=([^:]-[4])/a)o; [2] = [a/];y= anda) 0 isthe 

assumed frequency used in the calculation of the elements of 
the aerodynamic force matrix [A r ]. 

The solution of the above eigenvalue problem results in NM 
complex eigenvalues of the form 

<<V<o 0 K/y=p:±iv (7) 

The real part of the eigenvalue ( jl ) represents the damping 
ratio, and the imaginary part (v) represents the damped 
frequency; flutter occurs if jJ > 0 for any eigenvalue. 

In the present work, structural damping is set to zero and 
therefore damping is referred to as aerodynamic damping. 
Also, note that jl is opposite in sign to damping. 

TURBO Work-Per-Cycle Method 

The work-per-cycle approach is used to determine flutter 
stability. First, the flowfield through the blade row is 
calculated with no prescribed blade vibration. Starting with 
this converged steady flowfield, blade vibrations are 
prescribed in a selected mode, frequency, and nodal diameter 
pattern or phase angle. After the transients in the flowfield 
decay, and a periodic flowfield is obtained, the work done on 
the vibrating blade is calculated for a cycle of blade vibration 
as follows: 

W = § J- pdA-(8X/8t)dt (8) 

surface 

where, p is the blade surface pressure and A is the surface area 
vector. For harmonic vibration, the work-per-cycle of 
oscillation, can be rewritten as 
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j-pdA- 8q 0 G>cos(ei)t)dt (9) 

surface 

The aerodynamic damping ratio (Q associated with blade 
vibration is related to the work-per-cycle ( W) and the average 
kinetic energy {K F ) of the blade over one cycle of vibration 
through the following equation (Ref. 8) 



where 

K E =jj±mV 2 dt (11) 

t^C/Ccr ; C cr = 2/nco (12) 

In the preceding equations, C is the damping, C cr is the 
critical damping, m is the mass of the blade, V is the surface 
velocity due to blade vibration, and T is the time period. 

For small values of damping ratio which typically occur in 
aeroelastic calculations of interest, C,« 1, the aerodynamic 
damping ratio can be approximated as 

W/SkK e (13) 

If aerodynamic damping is negative, flutter can occur. Note 
that the structural damping (material and mechanical damping) 
has not been considered. Also, note that aerodynamic damping 
(Q is opposite in sign to jl . 


Results 

In this section, the results of the steady and unsteady 
computations are presented. Steady and unsteady 
computations were carried out at rotational speeds of 
100 percent (design), 85, and 75 percent. The part-speeds were 
selected to correspond to the speeds at which large vibratory 
responses due to flutter were encountered on the stall side of 
the operating line during testing. No flutter was encountered 
during testing at 100 percent speed. 

Note that in the present work, the harmonic balance 
aeroelastic code was run with very small amplitude vibrations 
of the fan blade to calculate linearized unsteady aerodynamics 
for a conventional linear flutter analysis. Also, the amplitude 
of vibration for the time-domain computations was in a linear 
range. 


Steady Computational Results 

In this study, the airfoil geometry used was calculated based 
on static deflections obtained from structural analysis for 
85 percent speed. The static deflections included the effects of 
rotational speed, applied pressures, and blade temperatures; 
nonlinear geometry effects were also included in the analyses. 
The blade geometry at 85 percent speed and nominal operating 
conditions was used for all computations. Previous 
calculations (Ref. 2) for this case have shown that changes 
from the nominal blade geometry due to changes in rotational 
speed were not significant and therefore a single geometry was 
used for computations at all speeds and operating conditions in 
the present work. The computational grid used was generated 
using commercial software. The grid for the harmonic balance 
computations is shown in Figure 1(a); the grid size is 
193x33x49 for the O-grid block that wraps around the blade 
airfoil with 193 grid points around the airfoil, 33 grid points in 
the circumferential direction, and 49 grid points in the 
spanwise direction. The H-grid blocks in the inlet and exit 
sections are each 17x33x49 with 17 points in the streamwise 
direction, 33 grid points in the circumferential direction, and 
49 grid points in the spanwise direction. The tip clearance was 
based on rig test measurements and is modeled using 9 points 
between the blade tip and casing. For the TURBO time- 
domain computations, the H-grid, shown in Figure 1(b), has a 
size of 121x51x39; 4 cells are used in the tip clearance region. 

The inlet flow conditions used in the computations 
consisted of circumferentially-averaged radial profiles of total 
pressure, total temperature, and flow angles. These prescribed 
profiles were based on rig measurements, supplemented by 
previous steady computations. The exit flow conditions 
consisted of a circumferentially-averaged radial profile of 
static pressure. This profile, which was based on pressure 
measurements at design speed, was used with uniform scaling 
for the computations at all speeds and at all conditions. 

Steady computations were carried out for three fixed 
rotational speeds (100, 85, and 75 percent). For each speed, 
the exit pressure profile was uniformly scaled to vary the mass 
flow rate and the operating point on the fan map. Thus, three 
speedlines were generated using each analysis code. Figure 2 
shows the convergence of the steady harmonic balance 
computations as a plot of non-dimensional torque with 
iteration counter (100 iterations per counter). The 
computations were carried out at 85 percent speed for various 
values of imposed backpressure at the exit boundary. As can 
be noted from Figure 2, excellent convergence was obtained 
for all operating conditions. Similar convergence was obtained 
for the computations at the other speeds as well as for the 
computations with the TURBO code at all speeds. All steady 
results presented here are from well-converged solutions. 
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Figure 1. — (a) Computational mesh for harmonic balance (HB) computations, 
(b) Computational mesh for time-domain TURBO computations. 



Iteration counter 


Figure 2. — Steady solution convergence (HB). 


The fan map is shown in Figure 3 with the results of 
harmonic balance steady computations denoted as HB and the 
results of the TURBO computations. Several computations 
were done for 75, 85, and 100 percent speeds. All computed 
results correlate well with measurements. Of particular interest 
is the last condition computed on the stall side, beyond which 
no converged solutions were obtained. For the HB results, the 
predicted pressure ratio is approximately 3 percent lower at 
75 and 85 percent speeds, and approximately 4.5 percent 
lower at 100 percent speed. For the TURBO results, the 



Figure 3. — Nondimensionalized fan map showing performance 
predictions from harmonic balance and TURBO codes. 

predicted pressure ratio is approximately 1.5 percent higher at 
100 percent speed and nearly identical at the other two speeds. 
For the HB results, the mass flow rate at the near-stall point is 
lower than measurements at 75 percent speed, in excellent 
agreement at the 85 percent speed, and significantly higher at 
the 100 percent speed. For the TURBO results, the mass flow 
rate is higher at the part speeds, and nearly the same at 
100 percent speed. Note that the stall line shown in Figure 3 
was estimated based on prior experience and was not 
confirmed during rig testing because of the occurrence of 
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flutter. Also, since the stall side was of primary interest, no 
attempt was made to determine the upper bound on choke 
flow. 

Overall, the HB results show good correlation with the 
TURBO results, except for a slight shift in the speed lines 
towards lower mass flow rate and lower pressure ratio for both 
75 and 85 percent speeds. At 100 percent speed, the HB 
results are shifted towards lower pressure ratio, and slightly 
higher mass flow rate. At 75 and 85 percent speeds, the HB 
calculations provide converged results for lower mass flow 
rates than were obtained with TURBO. Also, the HB results 
consistently show a slightly lower pressure ratio as compared 
to TURBO and the test data. These differences need to be 
investigated ftirther. It should be noted that significant 
differences exist between the numerical method used in the 
TURBO computations and the HB computations. These 
differences include time-domain versus frequency-domain, 
numerical discretization, grid, and algorithm used to solve the 
RANS equations, turbulence modeling, and others. Therefore, 
the differences in results noted in Figure 3 are not entirely 
surprising. 

To examine the sensitivity of the HB computed results to 
numerical parameters, the second and fourth order smoothing 
parameters, which are necessary to ensure numerical stability, 
were varied, respectively, between 0.25 and 1, and between 
0.002 and 0.006, as listed in Table 1. The results shown in 
Figure 4 indicate that no significant change in pressure ratio 
resulted from the stated variation in these numerical 
parameters. The mass flow variation was less than 0.5 percent. 


TABLE 1.— RANGE OF NUMERICAL SMOOTHING 
PARAMETERS USED TO STUDY VARIABILITY OF 
HB COMPUTED PERFORMANCE 


Case 

Smoothing parameter 


2nd order 

4th order 

1 

1 .0 (nominal) 

0.006 (nominal) 

2 

1.0 

0.002 to 0.005 

3 

0.25 

0.002 


1 


o 

-4— > 

03 


=3 

(/) 

(/) 


o 

o 

cr 


0 


Rotor map 
HB, 85% speed 



Inlet corrected flow 


Unsteady Computational Results 
Vibration Mode 1 

Unsteady computations were performed using the harmonic 
balance code at 100, 85, and 75 percent speeds to predict the 
aeroelastic stability of the blade. The first blade vibration 
mode and several nodal diameters covering the entire possible 
range were considered. Various steady conditions were 
considered from near peak-efficiency to near the stall line. 
This allowed the trend along the speed line to be assessed. 
Note that the results presented in this paper were all computed 
with only the 0th and 1st harmonics included in the 
aerodynamic analysis using the harmonic balance method. 

Figure 5 shows the convergence of the unsteady 
computations as a plot of non-dimensional generalized force 
(real and imaginary parts) with iteration counter (100 
iterations per counter) for various nodal diameters. These 
computations were for 85 percent speed and for an operating 
point near the peak-efficiency condition with a non- 
dimensional backpressure of 1.0. As can be noted, although 
the rates of convergence vary with nodal diameter, all the 
results are well-converged. Similar convergence 
characteristics were observed at the near-stall condition with a 
nondimensional backpressure value of 1.05. 

A sensitivity study was performed to investigate the effect 
of the unsteady grid scaling parameter. The nominal value of 
this numerical input parameter was increased by factors of 2 
and 5, and then decreased by the same factors. There was no 
change in the generalized force results within plotting 
accuracy - demonstrating that for the present calculations, the 
prescribed unsteady grid scaling has no effect on the results. 
Recall that in the present study, the harmonic balance 
aeroelastic code was run with very small amplitude vibrations 
of the fan blade to calculate linearized unsteady aerodynamics 
for a conventional linear flutter analysis. 


2xi o 4 



Figure 4.— Nondimensionalized fan map showing variability in Figure 5 ._hb unsteady solution convergence, 

results due to numerical smoothing parameters. 
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Figure 6. — HB generalized force variation with nodal diameter. 



Nondimensional backpressure 


Figure 7. — Variation of HB generalized force (imaginary part) 
with backpressure for A/D = 2. 




Nodal diameter Nodal diameter 


Figure 8. — (a) Variation of HB aerodynamic damping with nodal diameter at 85 percent speed, (b) Variation of TURBO 
aerodynamic damping with nodal diameter at 85 percent speed. 


Figure 6 shows the variation of the converged generalized 
force with nodal diameter for the two different values of 
backpressure at 85 percent speed. For a single vibration mode, 
the stability is determined by the imaginary part of the 
generalized force and Figure 6 shows that the stability varies 
significantly with nodal diameter. Further, with the change in 
backpressure from peak-efficiency towards stall, it can be seen 
that the imaginary part of the generalized force drops closer to 
zero at ND = 2. This trend is seen more clearly in Figure 7. 
The imaginary part of the generalized force is seen to clearly 
move towards zero as the backpressure increases - going from 
peak-efficiency towards stall. Note that Figure 7 includes 
results at nondimensional backpressure values of 1.052 and 


1.0522, with a change in the sign of the imaginary part 
occurring between these conditions. 

The generalized force was used to calculate the 
aerodynamic damping, which is plotted in Figure 8(a). The 
results are presented for backpressure values near peak- 
efficiency (1.0) and near-stall (1.05) conditions. The 
aerodynamic damping is seen to vary by an order of 
magnitude with nodal diameter of the traveling wave, with the 
minimum occurring at ND = 2 (forward traveling wave). Note 
that although the damping value is nearly zero at the near-stall 
condition, it is still positive even at its minimum value, 
indicating stability at a nondimensional backpressure of 1.05. 
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Nondimensional mass flow rate 

Figure 9. — Variation of HB aerodynamic damping with 
nondimensional mass flow rate for A/D = 2 at 85 percent 
speed. 

The aerodynamic damping calculated from the TURBO 
code is shown in Figure 8(b) as a variation with nodal 
diameter. Note that the backpressure values are slightly 
different for the TURBO results and the HB results. Also, the 
TURBO computations were done for fewer nodal diameter 
than with the HB code. A comparison of the results in 
Figure 8(a) and (b) clearly shows that both methods predict 
very similar variations with nodal diameter. Further the 
variation in damping values is nearly the same, as are the 
trends going from peak efficiency to near stall operating 
conditions (increasing backpressure). 

Additional calculations were carried out for a small increase 
in backpressure to identify the condition of zero aerodynamic 
damping (flutter). In the absence of these additional results at 
backpressure of 1.052 and 1.0522, it would be necessary to 
extrapolate the damping as a function of mass flow to identify 
the flutter point. Figure 9 shows the variation of aerodynamic 
damping with mass flow rate with a clear monotonic drop in 
stability along the speed line towards stall. 

For the TURBO computations, the aerodynamic damping 
was linearly extrapolated to obtain the flutter condition since 
the last computed condition near stall showed a small positive 
aerodynamic damping. 

Vibration Mode 2 

To verify that the mode 1 is the most unstable mode, HB 
unsteady computations were carried out for mode 2 at 
85 percent speed. The generalized force variation with nodal 
diameter, shown in Figure 10, is somewhat similar to that for 
mode 1 shown in Figure 6. Note that ND = 4 is the least 
damped nodal diameter. Unsteady computations were carried 
out for ND = 4 and 2 for various steady operating conditions 


4x10 4 



Figure 10. — HB generalized force variation with nodal diameter 
for mode 2 at 85 percent speed. 



Figure 1 1 . — Variation of HB aerodynamic damping with 
nondimensional mass flow rate for modes 1 and 2 at 
85 percent speed. 

(backpressure values or mass flow rate values). Figure 11 
shows that mode 1 2ND is the most unstable and mode 2 4ND 
has a very slightly higher damping. The rate of decrease of 
damping with decreasing mass flow rate is very similar for 
both modes as the stall line is approached. 

Other Rotational Speeds 

The results of unsteady calculations at 75 percent speed 
show trends that are nearly the same as those in Figures 5 to 9. 
Figure 12 shows the corresponding variation of aerodynamic 
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Figure 12. — Variation of aerodynamic damping with 
nondimensional mass flow rate for mode 1 A/D = 2 
at 75 percent speed. 



Nodal diameter 


Figure 13. — Variation of aerodynamic damping with nodal 
diameter for modes 1 and 2 at 100 percent speed. 



Nondimensional backpressure 


Figure 14. — Variation of aerodynamic damping with 
backpressure for modes 1 and 2 at 100 percent speed. 
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Figure 15. — Fan map showing calculated flutter point along 
with the measured flutter boundary. 


damping along the 75 percent speed line. Note that the 
aerodynamic damping is clearly negative at the lowest mass 
flow rate condition and the flutter point is obtained by 
interpolation. 

The results of unsteady calculations at 100 percent speed 
for modes 1 and 2 are shown in Figures 13 and 14. The 
variations with nodal diameter are somewhat similar to the 
those for 85 percent speed. However, the variations with 
backpressure (mass flow rate) show a very small slope in 
contrast to those seen in Figures 11 and 12. Based on the 
calculated aerodynamic damping and its small variation with 
backpressure, it is inferred that no flutter will occur at 
100 percent speed. 


The calculated flutter points at 85 and 75 percent speeds are 
plotted on the fan map along with the measured flutter 
boundary in Figure 15. It can be seen that the calculated flutter 
point is very close to the measured value at 85 percent sped 
with a difference of less than 0.5 percent in mass flow rate. As 
noted previously, the pressure ratio is slightly under-predicted 
and the difference from the measured data is 3 percent at the 
flutter point. For comparison, the flutter point calculated using 
TURBO is also plotted. Both calculations are very close to the 
test data; the present calculated result is slightly closer to the 
measurement in mass flow rate, but farther in pressure ratio as 
compared to the TURBO result. 
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At 75 percent speed, the difference between the present 
results and measurements is higher (about 4 percent) than at 
85 percent speed. The TURBO result also shows a larger 
difference from measurement at 75 percent speed than at 
85 percent speed. 

At 100 percent speed, no flutter point is predicted by the 
HB analysis and the flutter point predicted by TURBO is 
beyond the stall line, meaning that the operating condition will 
not be reached. 

Conclusions 

Steady and unsteady computations have been carried out 
using two Propulsion Aeroelasticity codes. These aeroelastic 
analysis codes solve the Reynolds-Averaged Navier-Stokes 
equations with blade vibrations using a time-marching method 
(TURBO) and a harmonic balance method (HB). The 
configuration selected was an experimental fan that 
encountered flutter during wind tunnel testing. The 
computational results were summarized and compared with 
experimental data. Overall, both TURBO and HB results are 
in good agreement with experimental data. The steady results 
were compared on the performance map and showed good 
correlation with data. However, the steady HB results under- 
predict the pressure ratio slightly and the steady TURBO 
results are closer to the measurements. Also, the HB analysis 
provides solutions at slightly lower mass flow rates than were 
obtained with TURBO at part-speed conditions. The HB 
results under-predict the pressure ratio by about 3%. A 
detailed look at the flowfields predicted by HB and TURBO 
may provide additional information regarding the source of 
the differences. A numerical study with the smoothing 
parameters in the HB analysis showed that the calculated 
performance is not sensitive to these numerical parameters. 

The flutter results from both HB and TURBO correlated 
very well with the experimental data. The least stable 
vibration mode matched experimental observations. The least 
stable nodal diameter pattern (interblade phase angle) also 
matched correctly with the experimental observations. The 
calculated flutter conditions agreed closely with the 
experimental data. A sensitivity study on the grid motion 
scaling parameter showed that the computational results were 
insensitive to variations in the selected value of this parameter, 
thus demonstrating the linearity of the results for the selected 
small amplitude of blade vibration. 

Computations for the first two vibration modes showed that 
the first mode had a slightly lower aerodynamic damping, 
correctly correlating with the experimental observations. 
Computations at three rotational speeds (100, 85, and 
75 percent) showed that flutter would be encountered at the 
part-speed conditions but not at the design condition, matching 
the experimental observations. 


In the present study, very small amplitude vibrations were 
prescribed to calculate linearized unsteady aerodynamics for a 
conventional linear flutter analysis as a first step to understand 
the characteristics of the harmonic balance code and of the 
flutter of an experimental fan. Future work will investigate the 
non-linear amplitude-dependent effects using the same 
harmonic balance aeroelastic code and the time-domain 
TURBO code. 
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